Electrostatic stability of insulating surfaces: Theory and applications 
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We analyze the electrostatic stability of insulating surfaces in the framework of the bulk modern 
theory of polarization. We show that heuristic arguments based on a fully ionic limit find formal 
justification at the microscopic level, even in solids where the bonding has a mixed ionic/covalent 
character. Based on these arguments, we propose simple criteria to construct arbitrary non-polar 
terminations of a given bulk crystal. We illustrate our ideas by performing model calculations of 
several LaAlC>3 and SrTiOs surfaces. We find, in the case of ideal LaAlC>3 surfaces, that cleavage 
along a higher-index (nlO) direction is energetically favorable compared to the polar (100) orienta- 
tion. In the presence of external adsorbates or defects the picture can change dramatically, as we 
demonstrate in the case of HbO/LaAlOa^OO). 

PACS numbers: 71.15.-m 



I. INTRODUCTION 

A uniformly charged plane separating two semi-infinite 
regions of space yields a divergent electrostatic energy; 
for this reason, polar surfaces or interfaces cannot exist. 1 
Yet, thanks to recent progress in epitaxial growth tech- 
niques, nominally polar terminations of insulating crys- 
tals are now routinely prepared and characterized within 
well-controlled experimental conditions'^ This is possi- 
ble because, in practice, there are several mechanisms 
available for a polar surface to neutralize the problem- 
atic excess charge, and possibly become thermodynam- 
ically stable. These include adsorption of foreign gas- 
phase species, changes in the surface stoichiometry, ionic 
and/or electronic reconstructions, or local metallization 
via accumulation of intrinsic free carriers; each of the 
above can, in principle, prevent the "polar catastrophe" 
by restoring the correct charge balance at the surface. 

Understanding and controlling these compensation 
mechanisms is a subject of great importance for many 
areas of fundamental science and technology. For ex- 
ample, surface polarity is of great interest for cataly- 
sis, gas sensing and energy applications^, as adsorp- 
tion and/or redox of gas-phase species are known to be 
strongly influenced by the electrostatic environment. In 
the context of perovskite-structure thin films and het- 
erostructures, control of surface charge/polarity is cur- 
rently investigated as a route towards the development 
of novel field-effect devices (e.g. in the LaA103/SrTi03 
system?—), or ferroelectric memories based on the tun- 
neling electroresistance effect^. Central to rationalizing 
all these phenomena is the intimate relationship between 
surface charge and bulk polarization in crystalline insu- 
lators, which was formally established by Vanderbilt and 
King-Smith in 1993»ii 

From the point of view of the theoretical analysis, it is 
crucial to establish an unambiguous criterion to classify 
a given surface as "polar-compensated" (i.e. an origi- 
nally polar surface that was neutralized via one of the 
aforementioned mechanisms) or intrinsically non-polar. 1 
This is not just a matter of nomenclature, but has very 



concrete practical relevance: non-polar terminations gen- 
erally tend to be more stable, as extrinsic (i.e. not orig- 
inated from the primitive building blocks of the insulat- 
ing bulk crystal) sources of compensating charge tend to 
have a high energy cost. Furthermore, if a given surface 
is polar, one needs to know precisely how much external 
charge is needed to neutralize it; this greatly facilitates 
the theoretical analysis by restricting the number of pos- 
sible candidate structures. We shall see in the following 
that, while the energetics is a genuine surface property, 
an exact answer to the latter question can be given al- 
ready at the bulk level. Many authors have already ad- 
dressed this issue in the past - we shall briefly mention 
hereafter the approaches that are most directly relevant 
to our work. 

Tasker 12 modeled a given ionic crystal as a lattice of 
point charges, corresponding to the nominal valence of 
the ions. Based on an abrupt truncation of this lattice, 
a given surface is then classified as polar on non-polar, 
depending on the behavior of the electrostatic energy. In 
particular, in the former (polar) case, the bulk repeated 
unit cell carries a finite dipole moment; this produces a 
diverging electrostatic potential unless compensated by 
an equal and opposite external surface charge density. 
This model, despite its simplicity, turned out to be sur- 
prisingly effective, and was able to correctly predict, at 
least at the qualitative level, the polar or non-polar na- 
ture of the vast majority of insulating surfaces. However, 
at the quantitative level this model has clear limitations. 
Many oxides and semiconductors display a marked cova- 
lent character, and the bulk polarization departs signif- 
icantly from the value that can be inferred from atomic 
positions and nominal valence charges. Hence the need 
for a more accurate treatment. 

To address these issues, and adopt a more realistic de- 
scription of the charge density of the solid, Goniakowski 
et aiX proposed a different criterion for classifying sur- 
faces as polar or non-polar. At the heart of the strat- 
egy of Ref. [l| is the concept of "dipole-free" unit cell. 
Given a certain plane orientation, one can demonstrate 
that it is always possible to choose a dipole-free repeated 
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unit along the normal to that plane; then, the remainder 
charge that is left at the surface (once the bulk units have 
been removed) determines the polar or non-polar charac- 
ter of the termination. This criterion, however, is not free 
from ambiguities: for the same surface termination there 
may exist more than one possible choice of the dipole-free 
bulk unit cell, which might lead to opposite conclusions 
about the polar or non-polar character of the surface (see 
section HTG2] for a detailed discussion). Also, identifying 
the dipole-free unit cell might be cumbersome in the case 
of higher index surfaces, where the structural complexity 
of the larger cell could complicate this type of analysis. 
Finally, the intuitive appeal of Tasker's model is appar- 
ently lost in the strategy of Ref. Q]: one needs to look at 
the ground-state charge density (e.g. as provided by a 
first-principles calculation) before drawing a conclusion. 

There are two further issues that are common to both 
methods. First, it is universally agreed that the surface 
polarity is a property of the actual lattice termination. 
This means that, for a given material and surface plane 
orientation, there might be polar and non-polar termi- 
nations, depending on the surface stoichiometry. How- 
ever, there is no established recipe to unambiguously de- 
cide, given a compound crystal and a surface orientation, 
whether a stoichiometric lxl non-polar termination is 
allowed at all. Furthermore, it is not clear how to con- 
struct, in general, a non-polar candidate structure with- 
out relying on a heuristic counting of the layer charges. 
Second, it was correctly recognized by both Tasker and 
Goniakowski that the issue of surface polarity is directly 
related to the bulk polarization of the material. How- 
ever, neither model traces a formal link to the modern 
theory of polarization in periodic insulators^, where the 
macroscopic P is a multivalued vector field, written in 
terms of the phases of the wavefunctions. Only the total 
charge density (modulus of the wavefunctions) is con- 
sidered in the model of Ref. [l], while explicit electronic 
orbitals are not addressed by Tasker's approach. Recent 
theoretical works have indeed highlighted the importance 
of the formal Berry-phase polarization in discussing po- 
larity at surfaces^ and interfaces^, but a general formu- 
lation of the problem, based on the formalism established 
in Ref. [ll], is still missing. 

Here we show that a Wannier function representa- 
tio n 15 i 16 together with the "interface theorem" of Ref. ITT1 . 
provide a very natural framework for addressing the 
above issues. Wannier functions were already shown to 
be a very useful tool, in layered superlattice o 17 ' 18 , for 
partitioning the polarization of a crystal into the con- 
tribution of individual charge-neutral units. Most im- 
portantly, Wannier functions are intimately linked to the 
modern theory of polarization in solids 19 , and therefore 
appear to be the most appropriate ingredient to discuss 
the issue of surface polarity, where the basic question con- 
cerns the existence of a finite dipolc moment perpendic- 
ular to the surface plane. We shall provide, based on this 
description, precise criteria to establish whether truncat- 
ing a bulk crystal along a given crystallographic orienta- 



tion can yield a non-polar surface. We shall demonstrate 
that answering this question involves only an analysis of 
the bulk, and that our scheme naturally leads to can- 
didate structures that can be used as a starting point 
for the subsequent determination of the thermodynamic 
ground state. To demonstrate our arguments, we focus 
on the surfaces of LaA103 and SrTiC>3, two prototypi- 
cal perovskite materials that have been at the center of 
the attention in the past few years as their polar (100) 
interface exhibits numerous peculiar properties. 

This work is organized as follows. In section |H] we 
introduce our definition of polar surface and its formal 
relationship to the theory of bulk polarization. We also 
establish a direct link to Tasker's model and we compare 
it to the "dipole-free" cell approach. In section UHl we ap- 
ply this formalism to a variety of systems, including non- 
polar LaAlG^nlO) and SrTiG^lll) surfaces. We also 
discuss electronic/ionic compensation mechanisms of po- 
lar LaAlO3(100). In section HVl we briefly address some 
related topics, including the case of ferroelectric surfaces, 
and possible extensions to covalent semiconductors. Fi- 
nally, in section [V] we present a brief summary and the 
conclusions. 



II. THEORY 

A. Definition of polar surface 

In full generality, for the surface of a crystalline insu- 
lator to be electrostatically stable, it must have a vanish- 
ing density of physical surface charge, tr sur f = 0. In order 
to introduce the notion of surface polarity, it is useful to 
separate er sur f into two distinct contributions, and rewrite 
the stability condition as 

CToxt + Pbulk ■ n = 0. (1) 

Here Pbulk is the bulk polarization, h is the normal to the 
surface plane, and cr cxt is a surface density of "external" 
compensating charges, which encompasses all contribu- 
tions that cannot conveniently be described as "bulk- 
like" in nature (we include the latter in Pbulk)- c cx t 
typically includes free charges (e.g. in the form of a con- 
fined electron gas) and/or bound charges (either in the 
form of surface adsorbates, vacancies, non-stoichiometric 
reconstructions, or non-isoelectronic substitutions). 

We define a given surface as non-polar if the stability 
criterion Eq. ([TJ can be satisfied in the absence of external 
charges er cx t, which implies 

Pbulk • n = 0. (2) 

This equation, at first sight, looks inconsistent with the 
current understanding of the surface polarity problem. 
It is now widely accepted that the polar or non-polar 
attribute is a property of the termination, not only of 
the material and surface plane orientation, contrary to 



3 



what Eq. ([2]) seems to suggest. We shall see in the fol- 
lowing section that the choice of the termination is only 
apparently absent from Eq. ([2]). It is implicitly included 
through the intrinsically multivalued nature of Pbuik) 
which is a well-established aspect of the modern theory 
of polarization in bulk insulators.— 



B. The bulk polarization 

1. As a Berry phase 

We consider a crystalline insulator described by three 
primitive translation vectors aj. 3 and a basis of N 
atoms located at positions R a , with a = 1,...,N. The 
"formal"— bulk polarization is usually defined as 

p ^^(E R ^- 2e E^fr)- ( 3 ) 

a— 1 i— 1 

Here Z a is the charge of the ionic core a, e is the (posi- 
tive) electron charge and (f>^ is the Berry phased along 
the reciprocal-space vector i; for simplicity we assume 
spin pairing, hence the factor of 2 in the electronic con- 
tribution. 

It is important to note that Pbuik i as 

defined in Eq. © 

is only defined modulo a "quantum of polarization" ; in 
other words, it is not a single- valued but a multi- valued 
function of the electronic and structural degrees of free- 
dom. This indeterminacy concerns both the ionic and the 
electronic parts in Eq. ©. On one hand, one has the free- 
dom to choose any of the periodically repeated images of 
each atomic specie, and thus change R Q by an arbitrary 
translation vector of the type AR = Uxaj + «2 a 2 + n 3 a 3- 
On the other hand, <f> \ are phases of complex numbers, 
and therefore only defined modulo 2n. 

In the following sections we shall use the equivalent 
formulation of Pbuik in terms of Wannier functions to 
illustrate the relationship between the multivaluedness 
of Pbuik and the termination of the crystal lattice. 



2. From the Wannier functions 

We shall explicitly assume, from now on, a single- 
particle picture, in terms of a Kohn-Sham set of orbitals, 
and a conventional (as opposed to topological) insulating 
state. Within these assumptions, it is possible to express 
the electronic ground state of the bulk solid in terms of 
a set of maximally-localized Wannier functions 15 , which 
are exponentially localized in direct spaced. Based on 
this representation, Eq. (J3J can be rewritten as 

1 N JV«i/2 

Pbuik- ^(E R «^- 2e E < r >;)' ( 4 ) 

a=l j=l 



where N c i is the total number of electrons in the primitive 
cell, and (r)j is the location of the center of the j-th 
Wannier function. 

Alternatively, we can think in terms of a charge density 
distribution that consists in the basis of atomic point 
charges together with the Wannier densities, 

N N.i/2 

Pcdi(r) = E Z <* 6 ( r " R «) ~ 2e E M r )| 2 ' ( 5 ) 

a=l j=l 

where Wj(r) is the J-th Wannier function of the prim- 
itive cell. By construction, the sum of all the periodic 
images of p ce ii "tiles" the total charge density of the ex- 
tended solid. Then, by combining Eq. (@| and Eq. ([5]) 
one immediately obtains the intuitive connection to the 
Clausius-Mossotti formula, 

Pbuik = ( 6 ) 

where d is the dipole moment of /9 C cii- 

This formulation provides a transparent way to par- 
tition the total charge density into individual primitive 
units, whose dipole moment correctly yields the formal 
value of Pbuik- In doing so, the phase indeterminacy of 
the electronic contribution to the polarization discussed 
in the previous Section has been reduced to a lattice in- 
determinacy, in all respects analogous to that character- 
izing the ionic contribution. In other words, all the com- 
plications related to the quantum-mechanical nature of 
the electrons have been mapped into a system of classi- 
cal point charges, where the atoms and the electrons are 
formally treated on the same footing. In the following 
section we discuss how this Wannier representation can 
be further partitioned into smaller units that retain the 
chemical information about the formal oxidation state of 
each ion, which is central to the notion of "polar surface" . 

C. Formal ionic charges 

The location of the Wannier functions generally re- 
flects the bonding properties of the material - in ionic 
solids they will cluster around the atoms, while in co- 
valent materials they will tend to occupy the bond cen- 
ters. We shall assume that the solid has at least a cer- 
tain degree of ionic character, so it is possible to "assign" 
each Wannier function to a given atom without ambigui- 
ties; this is certainly true in most known oxide materials. 
(With some caution the ideas developed here can be con- 
veniently adapted to any crystalline insulator; we shall 
briefly discuss the example of purely covalent semicon- 
ductors in Sec. IIVH We then combine each ion a with 
the Wannier orbitals j that "belong" to it and define a set 
of compound charge distributions that we call "Wannier 
ions" (WI), 

p§i(r) = Z a S(r) - 2e^ K(r + R Q )| 2 ). (7) 

jea 
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(We operated a translation so that the nucleus sits in 
the origin.) As the Wannier function locations usually 
agree remarkably well with chemical intuition, each of 
these N charge distributions will carry a monopole Q a 
corresponding to the "nominal" charge of the ion (e.g. 
™2e for O, +2e for Sr, +4e for Ti). In addition to their 
net charge, the WI are non-spherical and generally carry 
non-zero dipole moments d a . (Higher multipoles are also 
present, but are not directly relevant for the present dis- 
cussion.) p C oii(r) can now be rewritten in terms of the 
WI densities, 

p ccll (r) = 5>^(r-R Q ), (8) 

a 

which is equivalent to Eq. ([3]) except that here we use 
precautions to keep the basic WI units intact. It follows 
that the dipole moment of p cc ii(j) can be written in terms 
of two contributions, 

d = d PC +d W i. (9) 

The first term is the dipole moment of a system of point 
charges located at positions R Q , 

dpC — ^ I^ciQa- (10) 

a 

The second term is the sum of the individual dipole mo- 
ments of the WI, 

d wi = X! dct ' ( n ) 

Of- 
It is possible to show that d-wi is a single-valued, gauge- 
invariant quantity; this number contains all the non- 
trivial electronic contributions to the polarization that 
are due to the deformation of the ionic orbitals in the 
crystalline environment. The gauge invariance of dwi 
might be surprising at first sight, as the individual dipole 
moments d Q are manifestly gauge- dependent, i.e. they 
depend on the specific algorithm used to localize the 
Wannier functions. This arbitrariness cancels out when 
all d Q are summed up, as long as the assignment of each 
Wannier function to a specific lattice site remains unam- 
biguous. This is equivalent to stating that a given choice 
of R a uniquely determines the branch choice of the elec- 
tronic polarization, which is a reasonable assumption in 
ionic materials where the nature of the valence wavefunc- 
tion has typically a marked atomic character. 

With this new decomposition of d, we have overcome 
an important drawback of Eq. (0J and Eq. in the lat- 
ter two equations the decomposition of Pbuik into ionic 
and electronic contributions is physically meaningless - 
only the sum of the two terms is well defined. Here both 
dpc and dwi are formally meaningful objects. All the 
indeterminacy in the definition of d has been recast into 
the term dpc, which has an intuitive interpretation as 
the dipole moment of a system of classical ions, each of 
them carrying its formal valence charge. The term which 



contains the quantum-mechanical information about the 
electronic polarization effects is now dwij which is now a 
single- valued quantity; this has also a simple physical in- 
terpretation as the total dipole moment of the electronic 
clouds of the WI. 

Before closing this part, it is useful to point out a di- 
rect relationship between our formalism and the "layer 
polarizations" (LP), pi, introduced in Refs. [lj] and Il8l. 
If one is interested in layered perovskites with stacking 
axis along the (001) direction, it is useful to consider a 
decomposition of the charge density of the ABO3 cell into 
individual AO and BO2 layers. In the framework of the 
present work this implies grouping together the WI that 
belong to a given oxide layer. In particular, the total 
charge density of each layer I can be then written as a 
sum of all WI densities that belong to layer I, 

Mr)=5>$(r-R a ). (12) 

In II-IV perovskites the layers are formally charge- 
neutra l 17 ' 18 . pi(r) then carries a well-defined dipole mo- 
ment, which is directly related to the LP, 

Pi = J pi{z)zdz = y^(d a + RgQg) ■ z. (13) 

(The bar indicates in-plane averaging, and the integral 
is carried out along the stacking axis; S is the cell cross- 
section.) We shall illustrate this layer-by-layer decompo- 
sition of the total charge density with practical examples 
m section lHG~2l 

D. Crystal termination as a bulk property 

We shall illustrate in this section that the multival- 
uedness of the term dpc can be formally related to the 
surface termination of a semi-infinite crystal. This fact is 
not new, and was rigorously established within the mod- 
ern theory of polarization^. Here we discuss the impli- 
cations of the "interface theorem" for the electrostatics 
of polar surfaces. 

Following Goniakowski et aZ.— , we define a frozen bulk 
termination as a surface that is obtained by piling up 
"bulk unit cells" without any further electronic or ionic 
relaxation. For the time being, we shall limit our dis- 
cussion of surface polarity to this (somewhat unrealistic) 
type of surface, that we further specify hereafter; we shall 
make the link to more realistic surface models in the next 
section. In contrast with Goniakowski et a/.—, here we 
define our "bulk unit cell" as a charge density distribu- 
tion that results from a superposition of bulk WI, as in 
Eq. (JSJ. Then, we construct the charge density distribu- 
tion of the semi-infinite surface system as a superposition 

of /Ocell(r), 

p(r)= J2 Pcdi(r-R), (14) 

R-fi<0 
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where R = men + 1)23-2 +/?;$a3 is a real-space translation 
vector, and again n is the normal to the surface plane. 
This way of defining a frozen bulk termination has two 
crucial advantages: (i) The choice of using the compound 
WI as our "elementary particles" naturally ensures that 
every ion in the surface system will have exactly the same 
formal oxidation state as in the bulk. This is a central 
point in the definition of a polar surface - if we allowed 
for fractional orbital occupations no surface would be po- 
lar, (ii) The choice of cleaving the Bravais lattice, rather 
than the crystal lattice is particularly advantageous, as 
it naturally preserves the bulk stoichiometry everywhere 
in the system (if we allowed for stoichiometry changes of 
reconstructions again the notion of polar surface would 
be inconsistent). 

It can be easily verified that several types of unrecon- 
structed surfaces can be generated by using Eq. (TlH) . 
simply by changing the definition of p cc ii(r)- In partic- 
ular, we have the freedom to construct p C oii(r) in many 
different ways, simply by shifting each WI in the basis by 
an arbitrary Bravais lattice vector AR. Thus, the choice 
of the basis vectors r.; uniquely determines the surface 
structure, according to Eq. (|SJ) and Eq. (|T1| . (Of course, 
different choices of can lead to the same termination; 
the — ¥ termination" relationship is a many-to-one 
function.) On the other hand, we have shown in the pre- 
vious section that the choice of r; uniquely determines 
the value of Pbuik, out of the infinite possibilities allowed 
precisely by the arbitrariness in the choice of the basis 
vectors. This formally establishes the relationship be- 
tween Pbuik and the termination of the lattice. By con- 
struction, the physical net charge that lies at the surface 
of a frozen bulk termination as defined above is simply 
(T sur f = Pbuik • n, where Pbuik is the dipole moment (per 
unit volume) of an appropriate bulk unit cell, i.e. one 
that tiles the semi- infinite solid according to Eq. (fT4| . 

Within these assumptions, we define a frozen bulk ter- 
mination polar if the bulk building block used to con- 
struct it has a net dipole perpendicular to the surface 
plane; we define it non-polar otherwise. The problem of 
determining whether a surface is polar or not is, there- 
fore, reduced to the problem of calculating the dipole 
moment of a bulk unit cell made of WI. This, in turn, 
can be directly related to the result of a Berry-phase 
calculation in the bulk crystal, which can be routinely 
performed with most publicly available codes. In other 
words, the termination itself can be understood as a bulk 
property. 



E. Frozen and relaxed surfaces 

It might appear artificial to consider surfaces that are 
constructed by stacking electronic orbitals corresponding 
to bulk Wannier functions. At a real surface, electronic 
states always depart from their bulk counterparts be- 
cause of the peculiar chemical and electrostatic environ- 
ment produced by the truncation of the crystal. Further- 



more, also the ionic lattice undergoes nontrivial struc- 
tural relaxations in the surface layers, in response to the 
perturbation of the bonding network. A central point 
of our formalism is that both (electronic and ionic) sur- 
face relaxation effects are essentially irrelevant in the con- 
text of deciding whether a given surface is polar on non- 
polar. As a matter of fact, either type of relaxation only 
affects the surface dipole moment, and not the surface 
charge density. Thus, genuine surface properties (e.g. 
the alignment between the bulk bands and the vacuum 
levels, or the surface energy) certainly depend on these 
mechanisms, but the polarity (which depends only on the 
physical surface charge) won't be affected. This formally 
establishes the surface polarity as a property that can 
be completely understood at the bulk level - note that 
the termination dependence can also be understood as a 
bulk property as specified in the previous section. Then, 
all mechanisms that alter the surface charge (either in 
the form of a local composition change or as a modifi- 
cation of the formal oxidation state of the surface ions) 
are unambiguously understood as external compensation 
effects, and enter the definition of <7 ext . 

F. Construction of arbitrary non-polar 
terminations 

So far we have addressed the question of deciding 
whether a given surface, of a certain orientation and ter- 
mination, is polar or non-polar. One could wonder now, 
for a given bulk compound, (i) whether non-polar termi- 
nations can be constructed at all; (ii) if yes along which 
surface plane orientation; finally, it would be helpful to 
(iii) identify candidate non-polar surface structures based 
only on bulk information. In this section we shall illus- 
trate how this is done within the present definition of 
surface polarity. 

Essentially, the question (i) boils down to finding all 
possible values of Pbuik- This is, within the modern the- 
ory of polarization, a periodic lattice of points. The dif- 
ference between two arbitrary values of P is a multiple 
of a real-space primitive translation vector, 

p buik ^ p buik = ^(* a i +3&2 + fca 3 ). (15) 

Here Qq = ne is an integer n times the electron charge 
e. n, which determines the resolution of the Pbuik mesh 
depends on the convention of how the Wannier functions 
and the ion cores are grouped together. In particular, the 
constraint adopted here of assigning each Wannier func- 
tion to a specific ionic site generally restricts the lattice 
of possible values of Pbuik to a subset of those allowed by 
Eq. (@} and Eq. §5$ . This can be understood by observing 
that the new elementary building blocks of the lattice are 
the "compound objects" WI, rather than single electrons 
or ions. 

Answering question (ii) consists in finding the inter- 
sections between the infinite lattice of Pbuik values and a 
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given surface plane, which is a straightforward geometri- 
cal problem. 

Answering question (iii) then is easy by recalling the 
direct relationship between a given value of Pbuik and 
the dipole moment of a well defined bulk unit. More 
specifically, once a value (or a subset of values) of Pbuik 
is found for which Pbuik ■ n = 0, models of the non-polar 
surface can be readily built by stacking [using Eq. (1141) ] 
bulk unit cells that correspond to those same values of 
Pbuik- We shall present several practical examples of this 
strategy in Sec. Mil 



G. Relationship to previous approaches 

1. Tasker model 

Eq. ([5]) constitutes the rigorous link between Tasker's 
model' 1 - and the modern theory of polarization in peri- 
odic insulators. Within our formalism, the total excess 
charge at a frozen bulk termination can be exactly writ- 
ten as 



7 SU rf = Pbuik • n — 



(E a R aQa + d W i) • fl 



(16) 



where the sum is extended over all atoms in the semi- 
infinite crystal. The only difference between Tasker's 
model and Eq.[l6]is the additional, purely electronic con- 
tribution dwi, which comes from the polarization of the 
Wannier ions in the crystalline environment. This con- 
tribution vanishes in all solids that are characterized by 
a center of symmetry; in these materials the discussion of 
the surface polarity problem in terms of nominal charges 
is therefore rigorous and exact. Even in materials where 
dwi 7^ 0, neglecting this term is usually not crucial to as- 
sessing the polar or non-polar nature of a given surface. 
However, considering the WI contribution is essential for 
a quantitative estimation of a (which is the excess charge 
that needs to be compensated); this is especially true 
in ferroelectric materials, which generally have a large 
anomalous contribution to P. Thus, our formalism pro- 
vides a formal justification to Tasker's model, and com- 
pletes it by introducing an additional well-defined elec- 
tronic dipolar contribution, dwi- 

In addition to this, our strategy has important practi- 
cal advantages. Tasker's approach involves a direct cal- 
culation of Eq. (|16p by means of an infinite lattice sum, 
whose convergence is ensured by using Ewald summa- 
tion techniques. This procedure might be cumbersome 
in practice, and it requires a specialized computer code 
to perform the calculation. Our strategy greatly sim- 
plifies the problem, by reducing it to the calculation of 
the dipole moment of a small set of point charges. This 
can be done with paper and pencil in few minutes for a 
surface of arbitrary orientation, provided that one knows 
dwi- This vanishes in many cases of practical interest 
- whenever it doesn't vanish, only a single bulk Berry- 
phase calculation is needed to evaluate this contribution. 



Moreover, our strategy allows one to easily answer a num- 
ber of physical questions that were difficult to address 
within Tasker's approach, e.g. those discussed in the pre- 
vious section. 



2. Other approaches 

In order to fully appreciate the advantages of our for- 
malism, it is useful to compare it, in a practical case, 
with the alternative notion of "dipole-free unit cell" pro- 
posed by Goniakowski et alX For illustrative purposes, 
we consider the (100) surfaces of two prototypical per- 
ovskite materials, LaA103 and SrTiC>3, in their cubic 
high-symmetry phase. We shall use either the arguments 
developed in the previous sections, based solely on bulk 
properties of the respective materials, or the theory of 
Ref. llj, in order to assess the polar or non-polar charac- 
ter of these surfaces. 

In Fig. [1] we plot, along the (OOl)-oriented z axis, 
the xy-planar average of the total valence charge den- 
sity. (The left panels refer to LaA103, the right ones 
to SrTiC>3.) The upper panels show a possible decom- 
position of the electronic charge that leads to a dipole- 
free unit cell, which we construct as follows. First, we 
count the total charge of the ionic cores of the indi- 
vidual oxide layers. With the pseudopotentials used in 
this work, these are LaO(+17), A10 2 (+15), SrO(+16) 
and Ti02(+24). Next, we decompose the total valence 
charge by cutting it with abrupt (001) planes located in 
the interstitial regions. The location of those planes is 
chosen as to (i) respect the inversion symmetry of the 
crystal, and to (ii) assign to each layer an electron den- 
sity that exactly cancels the positive core charge of that 
layer. The resulting electron charge assigned to the AO 
layers is highlighted with a dark shading (light for the 
BO2 layers). By construction, the "unit cell" obtained by 
combining two adjacent layers (evidenced by the arrow 
and dashed lines in the figure) has zero dipole moment in 
both LaA103 and SrTi03. Hence, this construction fails 
at detecting any fundamental difference between LaA103 
and SrTiOs: both are predicted to have non-polar (001) 
surfaces. Of course, this prediction relies on a completely 
arbitrary partition of the total electronic charge density. 
There are many other ways to do it. For example, if one 
chooses a different location of the cut planes (e.g. at the 
mid-point distance between the atomic planes), or yet a 
more sophisticated prescription (e.g. based on the Bader 
analysis) , one generally gets a non- vanishing layer charge 
in both LaAlOs and SrTi03. From this perspective, one 
would have to conclude that the (001) surfaces of both 
materials are polar. The main point that we want to 
stress here is that, if we base our analysis solely on the 
total electronic density [as we have done in Fig. [Tfa-b)] , 
(i) the choice between one partitioning scheme and the 
other is arbitrary; (ii) any statement about the surface 
polarity inferred from such a partitioning is ambiguous; 
and (iii) such an analysis cannot be linked in any ways to 
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FIG. 1: (a-b) Decomposition of the total valence charge density according to the "dipole-free unit cell" picture. Light and 
dark shadings indicate the portions that belong to the BO2 and AO layers, respectively. We impose both distributions to be 
symmetric around the respective atomic layer locations, and to contain a number of electrons equal to the total valence charge 
of the ions. The unit cell indicated by the arrow and the dashed lines has zero dipole moment by construction in both LaA103 
(a) and SrTiC>3 (b). (c-d) Decomposition of the total valence charge based on maximally localized Wannier functions. The 
total Wannier densities of the AO and BO2 layers are shown by thick solid and thick dashed curves, respectively. The total 
charge density (obtained by the superposition of the periodically repeated Wannier densities) is shown as a thin solid line. 



the bulk polarization of the material (the latter cannot 
be defined, even in principle, in terms of the total charge 
density of a periodic crystal) . 

In Fig. QJc-d) we demonstrate how the Wannier-based 
decomposition of the valence density solves this problem. 
The thin solid lines show, as above, the ground-state elec- 
tronic charge densities p(z) (again, the bar symbol on 
p indicates that an in-plane averaging was performed). 
The total Wannier densities of each layer, defined as pi 
in Eq. (fT2|) , are shown as thick lines (solid for the AO lay- 
ers, dashed for the BO2 layers). Note that we show the 
electron density as positive, and we omit the bare pseu- 
dopotential charges from the plots. (These are a lattice 
of Dirac delta functions, centered at the oxide layer loca- 
tions.) As the Wannier functions are discrete objects, the 
total electronic charges are integer numbers. Most impor- 
tantly, the Wannier functions carry some crucial informa- 
tion (that is absent in the total valence density) on how 
the localized bound charges are organized in the insulat- 
ing state of each compound. It turns out that (summing 
up the contributions from the cores) the LaO and AIO2 
layers have a total charge of +1 and -1, respectively, while 
the SrO and TiC>2 layers result charge-neutral, in perfect 
agreement with the naive assumption of perfect ionic- 
ity Thus, the Wannier decomposition correctly identi- 
fies LaAlO3(001) as polar and SrTiC>3 as non-polar, in 
agreement with Tasker's classification. 



To corroborate our arguments, a further consideration 
is in order. One could be tempted to criticize our rea- 
soning by observing that Wannier functions are by no 
means uniquely defined starting from a given set of Bloch 
orbitals. In Fig. (Uc-d) we have chosen a maximally lo- 
calized^ representation, but there is nothing really fun- 
damental behind this choice. Is the identification of a 
surface as polar or non-polar robust against this arbi- 
trariness? The answer is yes. The formal proof of this 
statement was derived in 1993ii, several years before the 
maximally localized Wannier functions were first intro- 
duced. Our choice of a maximally localized representa- 
tion is motivated by its intuitive relationship to elemen- 
tary chemical concepts (e.g. formal valence), but the 
reader should keep in mind that this is just a convenient 
way of expressing a concept that has solid mathemati- 
cal grounds. In particular, the exact value of the surface 
charge can be computed at the bulk level, regardless of 
the degree of ionicity of the material. 11 
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(-1/2,-1/2) (-5/2,-1/2) (1/2,-1/2) (-3/2,-1/2) 




(0,0) (-2,0) (0,0) (0,0) 

FIG. 2: Relationship between Pbuik and the dipole moment 
of the bulk unit cell. Red balls represent O ions, with charge 
Qo = — 2e. Gold balls represent the A-site cation, either 
Sr (Qsr = +2e) or La (Qh a = +2e); green balls are the B- 
site cation, either Ti (Qti = +4e) or Al (Qai = +3e). The 
sketches (a-d) represent the projection of the atomic positions 
onto a (lOO)-oriented plane. On the top and the bottom are 
shown the values of Pbuik (in units of e/aj}, where ao is the 
lattice parameter in either material) resulting in LaAlOs and 
SrTiOa, respectively, from each cell arrangement. 



III. APPLICATION TO PEROVSKITE 
SURFACES 

A. LaAlOa and SrTiOa: bulk properties 

We now use two prototypical perovskite materials, 
LaAlOs and SrTiC>3, to illustrate our strategy in prac- 
tice. This choice of materials is motivated by the re- 
cent discovery of a conducting electron gas at their polar 
(100) interfaced This has generated a lively excitement 
in the research community and a renewed interest in the 
theoretical foundations of the surface/interface polarity 
problem^d 

We shall address questions (i-iii) raised in Sec. Ill F[ For 
simplicity, in both materials we consider only surfaces of 
the type (Oij). This means that only the projection of 
Pbuik on the yz plane is relevant, and our procedure can 
be conveniently represented on 2D graphs. Our strat- 
egy is general, and this choice was made only to simplify 
the notation and the graphical representations. We also 
consider both bulk compounds within their high symme- 
try cubic phase. (Both materials are characterized by 
zone-boundary distortions, related to rotations and tilts 
of the oxygen octahedral network; however, as these dis- 
tortions are non-polar in nature, they are irrelevant for 
the present discussion.) 

To start with (question i), we need to find the lattice 
of "allowed" values of Pbuik in either material. Recall- 
ing that both compounds are characterized by a center 
of symmetry, it follows that dvvi = 0, and Pbuik is ex- 
actly determined by the formal valence charges of the 
participating ions, all sitting in their high-symmetry lat- 
tice sites. Now, the formal ionic charges are La(+3), 
Al(+3) in LaA10 3 , Sr(+2), Ti(+4) in SrTi0 3 ; oxygens 
in either compound have a formal charge of (-2). In Fig. [2] 
we show how different choices of the crystal basis of five 
atoms lead to different dipole moments per unit cell, and 




FIG. 3: Lattices of allowed values of Pbuik in LaAlOa (left) 
and SrTiOa (right). In the left panel we show four possible 
surface plane orientations. The polar orientations (dashed red 
lines) have no intersections with the Pbuik lattice. The con- 
trary is true for the non-polar orientation (solid black lines). 



hence to a different formal polarization. If we could take 
all (infinite) combinations, we would obtain an infinite 
lattice of points, which is isomorphic with the real-space 
Bravais lattice of the cubic crystal. The 2D projection 
of the lattice of Pbuik in either compound is shown in 
Fig. [3] Even if the two compounds are isostructural (re- 
call that we consider both LAO and STO in their cubic 
phase), there are two important differences in their for- 
mal polarization lattice. First, Psto is centered in the 
origin, while Plao is centered in (1/2,1/2,1/2). Sec- 
ond, Psto has a coarser mesh than Plao _ the spacings 
are doubled because the constituent point charges are all 
even in the former. Note that both Pbuik lattices are 
centrosymmetric, consistent with the absence of a spon- 
taneous polarization in either material^ 

To answer question (ii), we need to find all possible in- 
tersections between a surface plane and the allowed val- 
ues of Pbuik; the projection of a few representative sur- 
face planes are plotted in Fig.[3^a). As it can be readily 
appreciated from the diagram, the aforementioned quali- 
tative differences between the respective Pbuik lattices of 
STO and LAO have important consequences on the elec- 
trostatics of the surfaces. In particular, in STO the origin 
belongs to the allowed values of Psto, an d any surface 
plane intersects the origin by construction; therefore, a 
non-polar surface of any possible orientation can be read- 
ily constructed [we shall illustrate the case of STO(lll) 
in Sec. IIII C] , Conversely, in LAO only specific plane ori- 
entations intersect the Plao lattice [note that the (100) 
orientation is correctly classified as polar] . In the follow- 
ing Section we shall consider a subset of these (infinite) 
possibilities, i.e. the vicinal (Oln) surfaces, where n is 
an arbitrary odd integer number. We shall focus on the 
lowest-index cases with n = 1,3,5. 



B. Vicinal LaAlO;? surfaces 

We first construct preliminary models for the (Oln) 
surfaces with n — 1,3,5. These are obtained by using 
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FIG. 4: Slab models for the vicinal LAO surfaces described 
in the text. Top: (Oil). Center: (013). Bottom: (015). Color 
code for atoms is the same as in Fig. [5] Thick dashed lines 
indicate the supercells used in the simulations. Thin lines are 
guides to the eye. Shaded areas highlight the basic primitive 
unit that was used to construct the slab model. Thick arrows 
indicate the dipole moment of the basic unit, parallel to the 
surface plane. 



Surface type 




Surface orientation 






(on) 


(013) 


(015) 


A 


1.93 


2.23 


2.15 


B 


1.93 


1.85 


1.83 


Cleavage 


3.86 


4.08 


3.98 



TABLE I: Calculated energy per area for the LAO surfaces 
described in the text. An ideal cleavage of the crystal is as- 
sumed to leave a pair of A and B surfaces. All values are in 
J/m 2 . 

slab geometries, within the supercell method. First we 
choose a unit cell with the appropriate translational peri- 
odicity in plane, and enough room along the out-of-plane 
direction to accommodate both the slab and a vacuum 
region (slab and vacuum thicknesses are treated as con- 
vergence parameters). Second, we tile the slab region 
with repeated copies of a well-defined primitive basis of 
atoms, which is chosen in a such a way that its dipole 
moment lies exactly parallel to the surface plane. (This 
implies that the choice of the basis depends on the sur- 
face orientation.) This procedure leads to the slab models 
sketched in Fig. @] 

The first observation is that all these surface models 
[except maybe the (011) case] present alternating LaO- 
type and A102-type terraces, and these terraces tend to 
grow wider and wider for increasing n. Note that [again, 
with the only exception of the (011) case], the construc- 
tion described above produces, in fact, two inequivalent 
surface structures for each orientation. In other words, 
the models of Fig.Uldo not enjoy inversion symmetry. We 
shall refer to these two surfaces as "type A" and "type 
B" , where type A presents LaO-type step edges and type 
B has AlO-type edges. Remarkably, it is easy to real- 
ize that one can change from A-type to B-type simply 
by displacing an oxygen atom from one step edge to the 
neighboring one. This way, starting from the "mixed" 
AB-type slabs Fig. @] one can readily construct pure AA 
or BB slabs. One can verify that the resulting AA and 
BB models do enjoy inversion symmetry. Since going 
from A to B preserves the bulk stoichiometry, this al- 
lows for a rigorous definition of the surface energy for all 
individual surface structures. 

In practice, in the simulations we use a slab thickness 
of approximately 4-5 LAO cells in each case, which is 
more than sufficient to obtain a well-converged value of 
the surface energy. The surface energy is defined as 

E sm l = 7^(£ slab - A^bulk), (17) 

where .Esiab and £buik are the relaxed total energy of the 
slab supercell and of LAO bulk, N is the total number 
of LAO units in the slab model, and S = a^\/\ + n 2 is 
the surface area in each case (ao is the lattice parame- 
ter of cubic LAO; the factor of 2 takes into account the 
fact that a slab has two surfaces). In Table Q] we report 
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FIG. 5: Total density of states of the LAO(013) slab models. 
Surface A (black curve), B (red dashed) and bulk (shaded 
area) are shown. 



the results. Comparing these values with previous liter- 
ature studies is difficult, as studies of vicinal perovskite 
surfaces are scarce. Only the lowest-index (Oil) surface 
type has been investigated to some extent, although we 
weren't able to find data specific to LAO. Concerning 
other perovskite materials, Eglitis and Vanderbilt 2 ^ re- 
ported an energy of 1.52 J/m 2 for an isostructural O- 
terminated SrTiO3(011) surface structure. The value we 
obtain for LAO, 1.93 J/m 2 , is somewhat larger but oth- 
erwise of comparable magnitude. Note that in the study 
of Ref. H3I a different (hybrid) functional was used - LDA 
might well overestimate surface energy values due to the 
well-known overbinding issues. 

It is interesting to note that in the case of B-type sur- 
faces the energy decreases slightly for increasing index n. 
We ascribe this behavior to the lower steps-to-terraces 
ratio in the (013) and (015) surfaces (undercoordinated 
step sites are likely to be less favorable). We consider 
it unlikely, however, that this energy be further reduced 
for n > 5. Increasing n would lead to larger and larger 
terraces that are locally charged [either of the LaO(+) 
or AIO2O type], and the electrostatic cost (roughly lin- 
ear in n) would eventually dominate over the step energy 
(proportional to 1 jn) in a way that bears many analogies 
to Kittel's theory of domain walls. Still, the increased 
stability of the vicinal (013) and (015) surfaces [compared 
to the (011) orientation] suggests that these geometries 
could be, in principle, fabricated under appropriate ex- 
perimental conditions. The simultaneous presence AO 
and BO2 domains appears promising for applications, e.g. 
in selective self-assembly of functional nanostructures, as 
it was recently shown in the case of SrTiOs- 2 ^ 

The above considerations on the energetics haven't an- 
swered an important question yet: how can we verify that 
these surfaces are indeed non-polar, consistent with our 
predictions? A useful indication comes from the density 
of states. If a surface is polar, then there is a need for 



compensation via additional charge carriers (either elec- 
tron or holes) that deplete or populate the energy bands 
of the crystal in proximity of the problematic termina- 
tion. This typically results in a metallic surface. Con- 
versely, if the surface is non-polar, the bulk-derived Wan- 
nier functions alone are sufficient to ensure electrostatic 
stability, and therefore the system can remain insulating. 
In Fig. [5] we show the total density of states extracted 
from a (013) (A- or B-type) slab model, compared with 
the bulk LAO density of states. In all cases there is a 
wide gap separating the unoccupied from the occupied 
states. This fact, together with the inversion symme- 
try and perfect bulk stoichiometry of the slabs, directly 
demonstrates that the surfaces are non-polar, and that 
every atom contributes with a total number of electrons 
that exactly corresponds to its formal ionic valence. Sim- 
ilar considerations apply to the (110) and (510) surface 
models (not shown). 

It is important to stress that, contrary to a common 
misconception, all these surfaces are perfectly stoichio- 
metric by construction, and they are non-reconstructed 
as they have the highest possible translational symme- 
try that is allowed by each plane orientation. It is often 
assumed that the only "legitimate" structures that can 
be named frozen bulk terminations are those that are 
obtained upon cleavage of the crystal lattice, i.e. pre- 
serving the integrity of the bulk-like atomic planes. This 
is, however, just a convention that has nothing funda- 
mental to it. We believe it is more practical to truncate 
the Bravais lattice instead. This automatically preserves 
stoichiometry and translational symmetry, and dramati- 
cally simplifies the description of surface electrostatics. 

As a final remark, it is fairly easy to realize that all the 
(Oln) surface models presented in this section are non- 
polar for any non-ferroelectric perovskite material (or 
for a ferroelectric one in its high-temperature symmet- 
ric phase). This can be simply understood by observing 
that, by replacing the cations in each bulk primitive basis 
(see Fig. [3J with those of a different charge family (i.e. 
I-V or II- IV), the dipole moment changes its magnitude 
but not its direction. Also, all the surface models enjoy 
inversion symmetry and have ideal bulk stoichiometry - 
the absence of a dipole moment normal to the surface 
plane is therefore automatically guaranteed. Therefore, 
many of the considerations made here in relationship to 
the specific LAO case are actually completely general, 
and apply to all cubic perovskite compounds. 



C. Non-polar SrTi0 3 (lll) surfaces 

To further illustrate our arguments, we move now to 
the case of SrTi03. According to our definition of polar 
surface, and by observing that the simplest choice of the 
SrTiOs bulk unit has zero dipole moment, one would 
conclude that in SrTiOs any surface orientation is "non- 
polar". To illustrate this point we shall consider here 
the (111) surface, which has been classified as polar by 
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FIG. 6: Primitive bulk SrTi0 3 unit cells for the model A 
(left) and B (right) (111) slabs. Sr atoms (large green circle) 
lie at the corners of the cube; O atoms (small red circles) lie at 
the face-centered sites; Ti atoms lie at the center of the cube. 
Both choices of the primitive unit have zero dipole moment 
along any direction. 



most authors. Note that this surface is indeed polar if 
one insists on terminating the crystal lattice with either 
a Ti or a Sr03 layer. Our prescription of cleaving the 
Bravais lattice and tiling it with well-defined bulk-like 
formula units is less restrictive, and allows for non-polar 
terminations as we shall see in the following. 

We build two inequivalent stoichiometric slab models 
(that we call A and B) by stacking the primitive build- 
ing blocks schematically shown in Fig. [5] It is easy to 
verify that both choices of the primitive cell have zero 
dipole moment. (Again, to compute the dipole moment 
we use the formal charges. This is substantiated by 
the Wannier-based decomposition described in section ITU 
which provides the formal link to the theory of polariza- 
tion in bulk solids M-) The primitive translation vectors 
of the supercell are (in units of the bulk equilibrium lat- 
tice parameter a = 7.275 a.u.) ai = (^/l/2, y/S/2, 0), 
a 2 = (yiA-VW^O) and a 3 = (0,0,10); the out-of- 
plane spacing a 3 was chosen in order to include a suffi- 
ciently thick vacuum region separating the repeated im- 
ages of the 10-layer slabs. As the slabs do not enjoy in- 
version symmetry (there are a total of four inequivalent 
surfaces in our simulations), we apply a dipole correc- 
tion in the vacuum layer to avoid unphysical macroscopic 
fields in the bulk region of the SrTiC>3 films. We use a 
regular (8 x 8 x 1) T-centered A-point mesh to sample the 
surface Brillouin zone, and we fully relax our structures 
within the symmetry constraints allowed by the surface 
composition. Note that the A-type slab preserves the 
point group of the bulk (111) orientation, while the B- 
type slab has a lower symmetry due to the presence of 
an incomplete oxygen plane on one side. 

In Fig. [7] we show the relaxed structures of the two slab 
models described above (the primitive unit of the super- 
cell was repeated three times in both in-plane directions 
to obtain a clearer view of the structure) . Henceforth we 
shall indicate Al, A2, Bl and B2 the four inequivalent 
surfaces, where A and B refer to the specific slab model, 
and 1/2 refer to the top/bottom surface, respectively. 
Al has a Ti03-type termination (i.e. an ideal Ti-type 
surface where the Sr atom has been removed from the 
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FIG. 7: Relaxed geometries of the two (111) slabs described 
in the text. The left structure corresponds to model A; the 
right one to model B. 



topmost Sr03 layer), and correspondingly A2 contains a 
Sr-type termination, where undercoordinated Sr atoms 
protrude from the underlying oxygen group. Bl has a 
supplementary O atom accommodated on top of an ideal 
Ti-terminatcd surface, and this atom forms a tetrahe- 
dron surrounding the topmost Ti atom. B2 has an O 
vacancy in the terminating SrC>2 layer; model B can be 
therefore obtained from model A simply by displacing a 
neutral SrO unit from the 2 (bottom) to the 1 (top) sur- 
face. Most of these surfaces were already considered in 
Ref.[26|, and indicated as "small unit-cell reconstructions" 
of SrTi03(lll). We note that surface reconstructions 
are typically associated with a reduction in the trans- 
lational symmetry group, which is not the case for any 
of these models. Therefore, we rather regard these as 
primitive, stoichiometric bulk terminations. Whatever is 
the nomenclature, the authors of Ref. [2(j| correctly rec- 
ognized the formal charge neutrality of these "valence- 
compensated" terminations. 

Similarly to the LaA103 case, we analyze the electronic 
properties of these surface models to verify their insulat- 
ing character. We plot in Fig.[8]the local density of states 
(LDOS) integrated on spheres of radius 3.0 bohr sur- 
rounding the Ti atoms. In the main panels we show the 
average Ti LDOS in the middle of the slab (gray shaded 
areas), which we take as our bulk-like SrTi0 3 reference 
curve. We also show the LDOS corresponding to the out- 
ermost Ti atom at the top (red dashed curve) and bottom 
(solid blue curve) surfaces. At Al the gap is smaller than 
in the bulk, as a narrow band of Ti-derived unoccupied 
orbitals splits from the conduction band. The band gap 
narrowing is rather extreme at A2, where a highly dis- 
persive surface state makes the gap as small as 0.1 eV 
at the r point (presumably this free-electron-like state 
is originated from the s and p state of the protruding 
Sr ions). To better illustrate this, we show a blowup of 
the LDOS in the inset. Here we also plot (thin black 
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FIG. 8: Local density of states (LDOS) on the Ti atoms for 
the two slab models described in the text. The Ti atom closest 
to top surface (type 1) corresponds to the red dashed curve; 
that lying closest to the bottom (type 2) surface is indicated 
as a solid blue curve. The average LDOS on the two central 
Ti atom of either slab is shown as a shaded gray area. Top 
panel: model A; the inset shows a blow-up of the spectrum in 
a neighborhood of the Fermi level (the LDOS of the Sr atom 
closest to the bottom A2 surface is also plotted as a thin solid 
curve). Bottom panel: model B; the inset shows the atomic- 
like spectrum of the topmost O atom (solid red curve with 
white shading), as compared with the bulk-like spectrum of 
an O atom lying far from the surfaces (dark areas) . Note that 
in the A case (top panel) we used a finer (16 x 16 x 1) fc-point 
grid to compute the LDOS, in order to better describe the 
dispersive surface state at the A2 termination. 



curve) the LDOS of the outermost Sr atom, where the 
surface state has its maximum weight. The nearly fiat 
DOS (the wiggles are caused by the finite k resolution) 
between and 1.5-2 eV, typical of a parabolic band in 
2D, is clear. The B slab presents, overall, an energy gap 
which is much closer to the bulk value; this suggests that 
the system is electronically more stable than in A. Both 
at Bl and B2 the gap reduction is caused by valence-band 
derived surface states; these are reminiscent of the states 
that are found at some B02-terminated (100) perovskite 



surfaces. Conversely, no conduction band-derived states 
are present. Especially interesting are the sharp peaks 
appearing at Bl; these are derived from the atomic-like 
orbitals of the outermost O atom. To illustrate this point, 
we plot in the inset the LDOS on the terminating O, 
which lies at the vertex of the surface tetrahedron sur- 
rounding Ti; for comparison, we also show the LDOS of 
a bulk-like oxygen far from the surfaces. The 2s- and the 
2p- derived features of the surface O appear extremely 
sharp and atomic-like, in contrast with the substantial 
broadening in the SrTi03 bulk caused by band disper- 
sion. A tctrahedral coordination might look unusual for 
Ti, which tends to adopt octahedral coordination in most 
(if not all) stable bulk oxide phases. Nevertheless, at 
the SrTiOa surface, analogous tetrahedral units were re- 
cently shown both experimentally and theoretically to be 
energetically favorable,— even in the case of the (110) ori- 
entation where there exist alternative (lxl) structures 
with relatively low energy.— 

Finally, we shall comment on the relative energy of 
these structures. Unlike the (nlO) models discussed in 
the previous section, here it is not possible to construct a 
stoichiometric and symmetric slab; therefore, we can only 
calculate a cleavage energy for A and B, E c \. Resolving 
this value into the contributions of the top and bottom 
terminations would require further considerations about 
the chemical potential of Sr, Ti and O; this goes be- 
yond the scopes of the present work. We find E C \(K)= 
6.27 J/m 2 and E cl (B)= 3.94 J/m 2 . These values are 
both larger than the previously reported cleavage ener- 
gies along the (100) or (110) directions. Especially the A 
model has a high energy cost, consistent with the "open" 
nature of the low-coordinated surface sites, and with the 
relatively unfavorable electronic configuration discussed 
in the previous section. Model B, on the other hand, has 
a cleavage energy that is significantly smaller, and (on 
average) reasonably close to typical (110) surface ener- 
gies. It should be kept in mind that the cleavage energy 
might be unequally distributed between the Bl and B2 
surfaces - we cannot exclude that one of the two might 
be quite stable in a wide range of thermodynamic con- 
ditions. -Ed(B) can be directly compared to the values 
reported in Ref. [26|, where it appears to correspond to 
the sum of the surface energies of model 3 and 4. The 
authors of Ref. [26| did not use LDA but a variety of dif- 
ferent density functionals, with values ranging from 4.94 
eV (PBE) to 6.41 eV (TPSSh) per surface cell. Our LDA 
value of 6.31 eV per surface cell compares favorably with 
the highest value reported there, consistent with the sys- 
tematic tendency of LDA towards overbinding. 

It is interesting to compare our calculated .E c i(A)=10.0 
eV/cell to the energy associated with the "textbook" 
cleavage, i.e. that leaving atomically flat, metallic and 
polar Ti / Sr03 terminations. Assuming that our LDA 
values are comparable to the TPSSh results of Marks et 
al., we can infer the Ti / Sr03 cleavage energy by sum- 
ming up the TPSSh surface energies of model 1 and 2 
in the aforementioned work; this yields a value of 12.3 
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FIG. 9: (a) and (c): Total DOS of the (100) LaAlOa slabs, 
highlighting the population of the states involved in the com- 
pensation of the surface polarity (shaded area, indicated with 
an arrow), (b) and (d): Plane-averaged density of compen- 
sating charge, related to the shaded portion of the DOS in (a) 
and (c). Top (a-b) and bottom (c-d) panels refer to the LaO- 
and A102-terminated slabs, respectively. 



eV / cell. As surprising as it may sound, our non-polar 
cleavage model A, with the severely undercoordinated Sr 
atoms protruding from surface A2, is still about 2 eV /cell 
lower in energy than the atomically flat cleavage model. 
This fact highlights the importance of achieving electro- 
static stability from bulk-like building blocks, without in- 
voking external compensation mechanisms (such as hole 
or electron doping as in the case of Ti / SrOs). 



D. Polarity compensation of LaAlOa(lOO) 

So far we restricted our analysis to ideal non-polar sur- 
faces, i.e. systems where only bulk-like building blocks 
are present. To complete our discussion we now con- 
sider a prototypical polar surface, LaAlO3(100), and il- 
lustrate how our arguments apply to the analysis of se- 
lected compensation mechanisms, where we introduce ex- 
trinsic sources of compensating charge <7 ex t • 



1. Via metallic earners 

First, we consider the clean (1 x 1) LaO- and AIO2- 
terminated (100) surfaces. Due to the built-in dipole of 
the bulk unit cell that we must use to construct these ter- 
minations, there is an excess charge of +0.5e and — 0.5e 
per surface unit cell, respectively. If we don't relax the 
translational symmetry, the only possible compensation 
comes from metallic carriers, either in the form of con- 
duction band electrons or valence holes. By performing 
two separate calculations of symmetrically terminated 



6.5-unit cell thick slabs, we indeed obtain metallic sur- 
faces. In Fig. [§] we plot the total density of states for 
both slabs, where the Fermi level clearly crosses either 
the valence band (AIO2 termination) or the conduction 
band (LaO termination). 

An interesting feature of the DOS of Fig.[9ja-c) is that, 
in both cases, a clear gap persists in the spectrum. This 
means that, in spite of the partial metallization, the con- 
duction and valence bands preserve their respective iden- 
tities. This observation implies that we can rigorously 
separate what we consider "bound charges" (which are 
all of bulk origin here, as we don't introduce extrinsic 
species in the supercell) from "external compensa ting 
charge" , following the prescriptions of Ref. [27| and [28|. 
The former, which we take as the total charge density 
of the (completely filled) valence-band manifold, are im- 
plicitly included in the definition of Pbuik! the latter can 
be either a positive external density of valence-band holes 
(ext,h) or a negative density of conduction-band electrons 
(ext,e), 

Pe*t,e(«) = / p(E,z)dE (18) 

* -Emid — gap 

/■-^mid-gap 

Poxt.hW = - / p{E,z)dE. (19) 

Here p(E, z) is the planar-averaged and energy-smeared 
local density of states defined in Ref. 127], and E-p is the 
Fermi level. The electronic states that contribute to the 
integrated charge densities p xt,c and p e xt,h are evidenced 
as shaded areas in the DOS plot of Fig. [SJ (Note that 
the DOS is the volume integral of p(E, z).) In Fig. [9jb- 
d) we plot the compensating surface densities p cx t,c and 
Pcxt.h- Both appear localized to the surface region, al- 
though they display a relatively slow decay into bulk 
LaA103, and amount (within machine precision) to a to- 
tal of exactly plus or minus half an electron per side. 
This demonstrates the full consistency [in the sense of 
Eq. ([I])] between the "external charge" defined in Eq. (JTSJ) 
and (UHJ), and the prediction of excess bound charge com- 
ing from the analysis of Pbuik- Note that, in the case 
of the LaO-terminated slab, part of the charge spills out 
into the vacuum region. This is a consequence of the vac- 
uum level being very close to the conduction band edge, 
which is populated by the compensating electrons. Con- 
versely, only 0(2p)-derived states contribute to p C xt,h- 
By combining the total energies of the reference slabs 
and subtracting an appropriate number of bulk reference 
units, we obtain a relaxed cleavage energy of 4.53 eV per 
surface cell (5.13 J/m 2 ). This is larger than the cleav- 
age energies we computed in Sec. IIII Bl for the primitive 
non-polar (nlO) surface models. 

2. Via external bound charges 

We shall now consider a different compensation mech- 
anism, where instead of metallic carriers the surface ac- 
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FIG. 10: Relaxed structure of the LaAlO 3 (100) slab with 
LaO (top) and AIO2 (bottom) terminations, compensated 
with OH(-) and H(+) groups respectively. 



quires bound charge via adsorption of external species. 
As a matter of fact, this is of concrete relevance for the 
interpretation of a vast number of physical phenomena. 
A real surface is always in contact with an atmosphere 
where various gas-phase species are present, and a num- 
ber of exchange/ adsorption/redox processes are usually 
thermodynamically accessible. On a more general basis, 
there are many situations where the surface layer dif- 
fers, either compositionally and chemically, from the bulk 
of the crystal. Consider, for instance, the Sr-decorated 
Si(100) surface that is used to promote coherent epitaxial 
growth of SrTi03i 29 ' 30 Given their technological and fun- 
damental importance, it is useful here to provide some 
examples, without the pretention to be exhaustive, of 
how our arguments can be translated to address those 
situations. 

In the specific case of LaA103, H atoms adsorbed at 
the surface of a film deposited on a SrTiO-3 substrate 
were found to significantly alter the electrical boundary 
conditions, e.g. by reducing or enhancing the residual in- 
ternal electric field in LaA103 and by influencing the free 
carrier concentration at the interfaced Experimentally, 
a humid atmosphere was demonstrated to be necessary 
to stabilize conducting paths at the buried interface^!. 
In turn, these "writing" and "erasing" processes^ appear 
to be mediated by charged surface adsorbates^. All in 
all, there is growing evidence that OH and H species are 
crucial to explain many outstanding phenomena exper- 
imentally observed at LaA103 surfaces and thin films. 
Hence the motivation for studying H20-based compensa- 
tion mechanisms, where the surface retains the insulating 
character of the bulk. 



Adding external species to a (1 x 1) surface won't make 
it insulating, as the excess charge is half an electron per 
cell (a periodic array of external species can provide only 
integer multiples of e). Doubling the surface unit cell 
leads to exactly plus or minus one electron, that now al- 
lows for an insulating state. We compensate this excess 
charge with a split water molecule (H adsorbed on the 
"negative" AIO2 side, OH on the "positive" LaO side) per 
\/2 x \/2 surface cell. We use a stoichiometric LaA103 
slab with a thickness of four unit cells, a c(2 x 2) in- 
plane translational symmetry and we relax the structure 
without imposing any symmetry constraint. As usual, 
we use a vacuum dipole correction to ensure the correct 
cancellation of the macroscopic electric fields due to the 
asymmetry of the slab. At equilibrium, the structure ap- 
pears as in Fig. [TO] Note the tilted position of the H 
atoms on the AIO2 side, consistent with the geometry 
found in Ref. for the H-LaA10 3 /SrTi0 3 system. The 
OH groups on the LaO side lie in a bridge site between 
two surface La atoms, thus occupying a natural lattice 
site for O. (An analogous location of the OH group was 
found on the SrO-terminated SrTi03 surface decorated 
with dissociated water—.) The system has a large in- 
sulating gap, almost equal to the bulk value, suggesting 
that this configuration might be fairly stable. We can 
estimate the energetics by considering a "wet cleavage" 
experiment where two LaA103 (001) surface are created 
and at the same time one free H2O molecule is split be- 
tween the two terminations, 



lab 



E 



H 2 0- 



(20) 



Here -E s i a b is the energy of the supercell described above 
with two LaA103 cells per surface unit and a thickness 
of four unit cells; -Ebuik is the bulk energy, calculated by 
including the antiferrodistortive tilt of the O octahedra, 
which now are allowed by symmetry; -Eh 2 is the energy 
of a free water molecule, calculated by using a cubic box 
of approximately 10 A lateral size. The resulting cleav- 
age energy per surface area is 2.50 J/m 2 , which is the 
lowest value calculated in this work. This result suggests 
that adsorption of OH and H groups, which are ubiqui- 
tous in most experimental setups, is a very likely candi- 
date to stabilize the LaA103 surface polarity. A study 
of LaAlO3(100) compensation via point defects was also 
recently reported in Ref. l33l . 

As a final remark, note that bound compensating 
charges, unlike the metallic carriers mediating electronic 
compensation, come in discrete units of e. Therefore, in 
cases where bound-charge compensation occur, it is most 
appropriate to "count" the external charges per unit area, 
which must satisfy the relationship 



bulk • n = - — 



(21) 



Here Q = ne (with n integer) is the formal oxidation 
state of the external defect or adsorbate, and S is the 
surface area per defect. Note that there exist defects 
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(e.g. transition metal cations) that are stable in several 
oxidation states; of course, the actual Q that occurs in 
the situation of interest must be used in Eq. (|2T|) . In 
case of doubt, the Wannier-based analysis of Sec. Ullcan 
be used to assess the formal oxidation state of a given 
defect. 



IV. DISCUSSION 

Here we shall put our results in the context of the 
current state-of-the-art in the field, especially regarding 
the fundamental theoretical understanding of the electro- 
static stability of insulator surfaces. 

A. Insulating nature of the surface 

As we mentioned several times when discussing our 
applications, it is likely that a non-polar (in the sense 
specified in this work) primitive surface will have a well- 
defined surface band gap. Here we shall further specify 
this point, to prevent dangerous generalizations. 

It is certainly true that a polar surface with all the 
atoms in their bulk oxidation state cannot exist. If we 
insist on keeping the local stoichiometry fixed, some of 
the atoms must change their valence in order to avoid a 
diverging electrostatic energy. In many cases this pro- 
duces partially filled electronic bands and a metallic sur- 
face. It is not difficult to imagine cases, however, where 
the surface atoms may change their oxidation state while 
preserving a gap in the spectrum. This would happen, 
for example, whenever the excess/defect charge amounts 
to an integer number of electrons, and there are ions in 
the lattice that have multiple stable oxidation states, e.g. 
most transition metals. Oxygen might also, in principle, 
change its formal valence from -2 to -1 to compensate a 
net surface charge - such a mechanism, stabilized via the 
formation of a peroxo bond, was reported in the case of 
a SrTi03 surface by Bottin et alS^. Therefore, a polar 
surface does not necessarily lead to a metallic surface. 
There are several other compensation mechanisms avail- 
able (often accompanied by a reduction in the transla- 
tional periodicity) that leave the surface insulating even 
without changes in the stoichiometry. We stress that in 
this latter case, however, the formal oxidation state of 
some atoms must change. 

Also the statement that non-polar surfaces are insu- 
lating is far from universal. Indeed, by replicating a suf- 
ficiently pathological choice for the bulk primitive unit, 
one might end up with a surface that has a very awkward 
bonding configuration. This might produce a dramatic 
departure from the bulk bonding environment, and in 
such cases it is well possible that one or more surface 
bands may close the band gap. An example of how this 
may happen is provided by our A2-type SrTiC>3 surface, 
where the band gap is reduced to a tiny value of about 0.1 
cV. Gap closure would also occur for our (nlO) LaA103 



surface models for a large enough n; at some point the 
electrostatic energy of the large LaO- and AlOvtype ter- 
races would become too large and eventually the gap 
would close. Note that surface relaxation usually helps 
stabilizing the truncated bonding network; in cases where 
our arguments would predict an insulating and non-polar 
surface, it is not infrequent to observe that a sizeable 
band gap opens only after full atomic relaxation. 

In conclusion, the relationship between electrostatic 
stability and insulating nature is certainly not a rigorous 
one. It is nonetheless a useful guideline, in the sense that 
if the surface bonding environment is not too pathological 
and the solid has a marked ionic character one usually ex- 
pects a non-polar surface to be insulating. Note that the 
insulating/metallic nature of a formally non-polar termi- 
nation is more an issue of chemistry than of electrostatics 
- it boils down to the chemical driving force for the atoms 
to preserve their bulk-like oxidation state. If the surface 
becomes metallic this will happen through a local rear- 
rangement of the electron cloud; the total surface charge 
won't change. 



B. Covalency arguments and "weak polarity" 

The (001) surfaces of II- IV perovskites such as SrTiC>3, 
BaTi03, etc. are classified as non-polar within our defini- 
tions, consistent with Tasker's assumption of formal ionic 
valence. Goniakowski, Finocchi and NogueraA challenged 
Tasker's classification by invoking covalent bonding ef- 
fects, which would produce a smearing of the electron 
cloud. This, in turn, would produce a distribution of the 
charge between the oxide layers that differs from their 
formal ionic charges. According to this interpretation, 
SrTiO3(001) was classified as "weakly polar". 

The disagreement between the two interpretations is 
rooted in the way the electronic charge density is parti- 
tioned into individual building blocks at the bulk level. 
Choosing Mulliken or Bader populations inevitably leads 
to a fractional charge per ion that is typically smaller 
than the (integer) formal oxidation state, and the indi- 
vidual SrO and Ti02 layers appear charged. While these 
ideas have certainly some merit, there are severe draw- 
backs as well. The most important one is that, by insist- 
ing on a partition of the electron cloud based exclusively 
on the total charge density p(r), one thwarts any further 
attempt at linking the discussion of surface electrostat- 
ics to the theory of polarization in bulk solids. Powerful 
and rigorous results of the modern theory of polarization, 
e.g. the "interface theorem" , are inconsistent with a de- 
scription of surface polarity in terms of Mulliken, Bader 
or even Born effective charges. (The fact that in II-IV 
perovskites the individual AO and BO2 layer do not sat- 
isfy the acoustic sume rule separately is sometimes taken 
as a further argument in support of the weak polarity 
concept.) The theory of bulk polarization implies a wave- 
function-based partition of p{r). The natural tool in that 
sense are the spatially localized Wannier function as we 
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discussed in section [II] Covalency effects are irrelevant 
in this context, except that they are implicitly accounted 
for through the location of the ground-state Wannier cen- 
ters and their spatial spread. Using Wannier functions 
might appear unnatural and complicated at first sight, 
but eventually they really lead to drastic simplifications 
and to an intuitive physical picture. In particular, if one 
wants to recover the intuitive classical formula P • n = a, 
there is simply no other way around. 

We believe that full consistency between the theory of 
bulk polarization and the theory of surface polarity is a 
must. Therefore, we caution against the use of concepts 
such as weak polarity or covalent charges as they are are 
inconsistent with the former. 



C. Other oxide surfaces 

1. Ferroelectric perovskites 

As the concept of surface polarity is intimately linked 
to the polarization of the bulk solid, Pbuik, it is par- 
ticularly insightful to discuss cases where Pbuik has a 
non-trivial behavior, as in ferroelectric perovskite mate- 
rials. Consider a (lOO)-oriented slab of BaTi03, with the 
spontaneous polarization vector, Ps oriented along the 
normal to the surface. Imagine that we have a stoichio- 
metric slab with ideal BaO and Ti02 terminations, and a 
monodomain state with perfect lxl periodicity; assume 
also that P points towards the TiCVtype surface. 

Both surfaces are polar, as Pbuik ■ fi — Ps 7^ 0. Con- 
trary to the LaA103 example, however, here Pbuik • n is 
not a simple fraction (plus or minus one half) of the polar- 
ization quantum e/S. Here Ps = pe/S, where p is a real 
number of the order of 0.25-0.35 (depending on the in- 
plane strain imposed to the film). Therefore, it might be 
technically difficult in a calculation to construct a com- 
mensurate supercell where Ps is accurately compensated 
by an appropriate coverage of charged adsorbates or de- 
fects (unless p happens by accident to be exactly equal 
to a rational number with a small denominator). A pos- 
sible trick to circumvent this difficulty is using the so- 
called "virtual crystal approximation" (VCA) 2 -. Here a 
fractionally charged pseudopotential is introduced at the 
surface to reproduce the effect of a disordered array of 
defects with the appropriate coverage. This way, the sur- 
face can be made insulating and charge-neutral at a low 
computational cost; the price to pay is that the VCA does 
not lend itself easily to the calculation of surface-specific 
properties, e.g. the energetics of a given compensation 
mechanism. 

A second important example is that of I-V ferroelectric 
perovskites, e.g. KNb03. Here the ferroelectric contri- 
bution to the polarization, Ps, adds up to the "composi- 
tional" built-in cell dipole,— which is Pq = ±e/2S (as in 
LaA103, the layers are formally charged, although here 
AO layers are negative and BO2 are positive). Note that, 
if |Pg| were (again, by accident) equal to half a quantum 



of polarization, both Nb02- and KO-terminated (100) 
surfaces would be non-polar (i.e. Ps would cancel out 
Po), provided that the spontaneous polarization points 
in the correct direction. This would be away from the 
surface for the p-type Nb02 termination, and towards 
the surface for the KO termination. 

Finally, it is worthwile mentioning the case of BiFe03. 
This material appears complicated at first sight, because 
of the tilted polarization axis [Ps is oriented along the 
(111) direction] and the compositional layer charges of 
±1 (both Fe and Bi are formally 3+ ions). However, 
within the formalism established in this work, predicting 
the excess surface charge density that will be present at 
a given ideal termination becomes trivially simple. For 
example, at the FeC>2 (100) termination we have [just like 
in the case of the A102-terminated LaAlO3(100) surface] 
an excess built-in charge of —e/2S ~ —0.5 C/m 2 . To 
this value we need to add the projection of Ps along the 
(100) axis, which amounts to approximately the same 
value.' 55 Therefore, if Ps points towards the FcOo-type 
(100) termination, this surface will be in practice only 
very weakly charged or even neutral. It goes without 
saying that the composition of the surface "pins" the out- 
of-plane component of the polarization to a fixed value 
that cannot be switched (unless the surface composition 
itself is changed, see Ref. \3w . 



2. ZnO 

The use of ZnO in many technological areas, as well as 
the recent progress in fabricating tailored nanostructures 
and functional surfaces with this material, have gener- 
ated a widespread interest in the fundamental proper- 
ties of its polar (0001) surfacei 3 ' 37 i 38 . Several possible 
compensation mechanisms involving, e.g. metallic free 
carriers^!, hydroxy lation/protonation 3 ^ or stoichiometry 
changes 3 - have been proposed over the years. In spite 
of this activity, the question of exactly how much excess 
charge is present at the polar Zn- or O- terminated sur- 
faces is still a source of confusion. 

For instance, there is a common belief that, starting 
from an ideal unreconstructed termination, removal of 
1/4 of the surface ions will lead to perfect compensation 
of the polarity^ This would be true if ZnO crystallized 
in zincblende phase. However, bulk ZnO is wurtzite- 
type, which means that on top of the compositional 
(zincblende-like) dipole it has also a non-trivial sponta- 
neous Ps^2 This Ps is of course not switchable, unlike the 
ferroelectric materials discussed in the previous section, 
but it does need to be taken into account when com- 
puting the surface charge. First-principles calculations 
of Ps have reported relatively small values (compared to 
a hypothetical zincblende reference structure) of Ps in 
bulk ZnO, of the order of 0.02-0.07 C/m 239 ' 40 . This im- 
plies that the necessary correction to the zincblende-like 
excess charge of 0.5e per surface cell is of the order of 
0.01-0.03 electrons. Even if this correction is not large, 
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one should keep in mind that, in a hypothetical free-slab 
calculation of ZnO where 1/4 of the O (and Zn) surface 
ions have been removed, after full relaxation there will be 
a non-zero residual macroscopic electric field in the slab 
of approximately £ s i a b = Ps/(eoe T ). Here eo is the vac- 
uum permittivity and e r is the static dielectric constant of 
ZnO (including piezoelectric effects). In fact, this obser- 
vation was used to calculate the spontaneous polarization 
of wurtzite BeO several years before the modern theory 
of polarization was developed4i 



D. Semiconductor surfaces 

While our arguments apply most naturally to ionic 
materials, where the assignment of the localized Wan- 
nier charges to a given atom is unambiguous, with some 
care they can be easily adapted to covalently bonded in- 
sulators. The main difficulty is that in semiconductors 
(e.g. Si) the maximally-localized Wannicr functions tend 
to occupy bond-centered sites, and are shared between 
two atoms - assigning a given Wannier function to ei- 
ther atom that participate to the bond is then entirely 
arbitrary. Nevertheless, one can usually establish a rea- 
sonable convention for partitioning the bulk solid into 
well-defined units. For instance, in Si one could assign 
four spin-up Wannier functions (their centers would form 
a tetrahedron around the nucleus) to one atom and four 
spin-down Wannier functions to the other atom in the 
basis. (It might appear somewhat artificial to use such a 
spin-split basis; however, for the present discussion, the 
information about the spin is irrelevant, only the charge 
density of the Wannier functions really matters.) Then, 
this decomposition yields a basis of two Wis that individ- 
ually retain the full symmetry of the lattice, are charge- 
neutral and have zero dipole moment. Primitive Si sur- 
faces are then predicted to be non-polar, but chemically 
they will be highly reactive because of the singly occu- 
pied "dangling bonds" ; this picture is consistent with the 
widely accepted understanding of Si surfaces. It is easy 
to see that by saturating these bonds with H one always 
obtains a non-polar and chemically stable surface (H does 
not add a net charge density as it contributes one electron 
and one proton to each dangling orbital). Alternatively, 
one could supply one Sr atom every two dangling bonds; 
this stabilization mechanism is important for growth of 
perovskite oxide films on Si substrates i 29 i 30 Interestingly, 
in the case of the Sr-decorated surface, further oxidation 
does not change the surface charge count j2H as additional 
O atoms achieve a closed-shell configuration by incorpo- 
rating the electron pairs already present in the saturated 
dangling bonds. This is a system where oxygen adsorp- 
tion does not change the surface charge, in striking con- 
trast with typical ferroelectric surfaces^ 

Of course, one could prefer to use other conventions, 
e.g. assign two doubly occupied Wannier functions to 
each Si atom. This way the Si(100) or (111) surfaces 
would be understood as "polar" , and they are indeed 



polar if one insists on counting electrons two-by-two (the 
dangling bonds would need to be either empty or sat- 
urated, without the necessary countercharge to balance 
the electrostatics). This means that the concept of po- 
lar surface becomes somewhat ill-defined if the solid has 
no ionic character whatsoever. Note that the formalism 
developed in Ref. [Til on which the present work heav- 
ily relies, provides always a rigorous means of calculating 
the surface charge from bulk properties, regardless of the 
(ionic or non-ionic) nature of the insulator, and indepen- 
dently of the convention that one uses to "assign" the 
bound electron charges to a given lattice site. A more 
extensive treatment of the covalent case can be found in 
Ref. [ll|, and was recently discussed also in Ref. I22I 



E. Interfaces 

In this work we decided to focus on surfaces, which 
is a special case of interface between two materials (one 
of them is vacuum). Whenever the second material is 
another crystalline insulator, the same arguments apply, 
but the "electrostatic phase diagram" can be substan- 
tially richer. The simplest case is that of two materials 
that have the same crystal structure, and we assume co- 
herent epitaxy, i.e. both semi-infinite regions have the 
same in-plane periodicity and the same crystallographic 
orientation of the atomic planes. However, there might 
be more complex cases - for example, the participating 
materials have different bulk structures, or they are not 
oriented along the same crystallographic direction. In 
any case, the electrostatics is always governed by the in- 
tuitive classical formula, 

(Pa - Pi) • h = cText- (22) 

Here Pi, 2 is the polarization in either material, calcu- 
lated by choosing a certain basis for the primitive basis 
of atoms and Wannier functions; <7 C xt is the "remainder" 
interface charge, that is left behind once one removes all 
the bulk-like primitive units on either side; h is the nor- 
mal to the surface plane. As in the case of surfaces, we de- 
fine an interface non-polar if, for an ideal termination of 
both materials with the maximum allowed translational 
symmetry one has <7 ex t = 0. 



V. CONCLUSIONS 

In summary, we have revisited the concept of polar 
surface within the context of the modern theory of bulk 
polarization. Our definition, which is consistent with 
the bound (and discrete) nature of electrons in the in- 
sulating state of matter, puts Tasker's classification on a 
firmer theoretical grounds, and corroborates it at the mi- 
croscopic level. We further complete Tasker's formalism 
with an additional term, which comes from the polariza- 
tion of the electron cloud in solids that spontaneously 
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break space inversion symmetry. Our calculations of 
non-polar LaAlO^nlO) and SrTiC>3(lll) surfaces, and of 
compensation mechanisms at LaAlO3(100), demonstrate 
that our formalism provides a convenient way of describ- 
ing the net surface charge in terms of bulk polarization 
and external sources (either "bound" or "free"). We have 
also illustrated some practical analysis tools that can be 
used to monitor the equilibrium distribution of compen- 
sating charge in a calculation. We hope that these tech- 
niques will be helpful for future first-principles studies, 
and more generally as a conceptual basis to rationalize 
the many interesting phenomena occurring at the sur- 
faces of insulating materials. 
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